Molecular hydrogen in silicon: A path-integral simulation 
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Molecular hydrogen in silicon has been studied by path-integral molecular dynamics simulations 
in the canonical ensemble. Finite-temperature properties of these point defects were analyzed in the 
range from 300 to 900 K. Interatomic interactions were modeled by a tight-binding potential fitted 
to density-functional calculations. The most stable position for these impurities is found at the 
interstitial T site, with the hydrogen molecule rotating freely in the Si cage. Vibrational frequencies 
have been obtained from a linear-response approach, based on correlations of atom displacements 
at finite temperatures. The results show a large anharmonic effect in the stretching vibration, w 8 , 
which is softened with respect to a harmonic approximation by about 300 cm -1 . The coupling 
between rotation and vibration causes an important decrease in w s for rising temperature. 

PACS numbers: 61.72.jj, 61.72.uf, 63.20.Pw, 71.15.Pd 



I. INTRODUCTION 

Hydrogen can be incorporated into semiconductors 
both intentionally and unintentionally during manufac- 
turing processes carried out for technological applica- 
tions. It appears in these solids in a number of differ- 
ent configurations: as an isolated interstitial, bound to 
impurities, bound to native defects, in molecular form, 
etCfi^ In the early 1980s, isolated hydrogen molecules 
were predicted to be stable in crystalline semiconductors 
and to play an important role in the diffusion of hydrogen 
in these materials 4£ However, they were not unambigu- 
ously detected by spectroscopic methods until more than 
ten years later.^^ 

Vibrational transitions have been reported for intersti- 
tial H2 in Sir^ Ge^ and GaAs£ In these semiconduc- 
tors, theory predicts that the H2 molecule is stable at an 
interstitial tetrahedral (T) site and behaves as a nearly 
free rotatori 10 ' 11 This gives rise at low temperatures to 
two stretching local vibrational modes originating from 
para and ortho nuclear states, which are split due to ro- 
vibrational coupling^ 

Here we will concentrate on isolated hydrogen 
molecules in the bulk of crystalline silicon. The interest 
of this problem is twofold. On one side, it is important as 
a point defect in semiconductor physics, for its relevance 
in the hydrogen diffusion and stability in these materials. 
On the other side, from a fundamental point of view, H2 
in silicon is an example of a light molecule sitting and 
moving in a confined geometry, and one can study its be- 
havior when localized in a spatial region with extension 
of a few A. 

Earlier theoretical studies of molecular hydrogen 
in semiconductors have concentrated on determin- 
ing the lowest-energy site and stretching frequency 
of the molecule, including in some cases anhar- 
monic effects derived from the calculated potential- 
energy surface ; 10 ' 11 ! 13 ! 14 ! 15 as well as the quantum rota- 
tion of H2 molecules i 16 ' 17 Density-functional electronic- 
structure calculations in condensed matter are very reli- 



able, but they treat atomic nuclei as classical particles, 
and typical quantum effects like zero-point vibrations are 
not directly accessible. These effects can be included by 
employing harmonic or quasiharmonic approximations, 
but are difficult to take into account when large anhar- 
monicities are present, as can happen for light impurities 
like hydrogen. 

To consider the quantum character of the nuclei, the 
path-integral molecular dynamics (or Monte Carlo) ap- 
proach has proved to be very useful. A remarkable ad- 
vantage of this method is that all nuclear degrees of free- 
dom can be quantized in an efficient manner, thus in- 
cluding both quantum and thermal fluctuations in many- 
body systems at finite temperatures. In this way, Monte 
Carlo or molecular dynamics sampling applied to evalu- 
ate finite-temperature path integrals allows one to carry 
out quantitative and nonperturbative studies of highly- 
anharmonic effects in solids . 18 ! 19 

In this paper, the path-integral molecular dynamics 
(PIMD) method is used to study interstitial hydrogen 
molecules in silicon. Special attention has been paid to 
the vibrational properties of these impurities, by con- 
sidering anharmonic effects on their quantum dynam- 
ics and the ro-vibrational coupling at different tempera- 
tures. The results of the present calculations show that 
anharmonic effects lead to a significant decrease of the 
vibrational frequencies of the impurities, as compared to 
a harmonic approximation. We have analyzed the iso- 
topic effect on structural and vibrational properties of 
these molecules, by considering also molecular deuterium 
(D2). Path-integral methods analogous to that employed 
in this work have been applied earlier to study hydrogen 
in metals^ and semiconductors , 20 i 21 ' 22 i 23 ' 24 as well as on 
surfaces . 25 ' 26 In connection with the behavior of molecu- 
lar hydrogen in confined regions, H2 has been studied in- 
side carbon nanotubes by diffusion Monte Carlo^ Also, 
path-integral simulation methods have been extensively 
applied to study condensed phases of hydrogen in molec- 
ular form^^l 

The paper is organized as follows. In Sec. II, we de- 
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scribe the computational method and the models em- 
ployed in our calculations. Our results are presented in 
Sec. Ill, dealing with the kinetic energy of the molecules, 
spatial derealization, interatomic distance, and vibra- 
tional frequencies. Sec. IV includes a discussion of the 
results and a summary. 



II. COMPUTATIONAL METHOD 
A. Path-integral molecular dynamics 

In the path-integral formulation of statistical mechan- 
ics employed here, the partition function is evaluated 
by a discretization of the density matrix along cyclic 
paths, consisting of a finite number P (Trotter number) 
of "imaginary-time" stepSi 32 > 33 In the implementation in 
numerical simulations, this discretization gives rise to 
the appearance of P "beads" for each quantum parti- 
cle. These beads can be formally treated as classical 
particles, so that the partition function of the original 
quantum system is isomorph to that of a classical one. 
This isomorphism is obtained by replacing each quantum 
particle by a ring polymer consisting of P classical parti- 
cles, connected by harmonic springs, In many-body 
problems, the configuration space is usually sampled by 
Monte Carlo or molecular dynamics techniques. Here, we 
have employed the PIMD method, which has been found 
to need less computer time for the present problem. We 
have used effective algorithms for performing PIMD sim- 
ulations in the canonical NVT ensemble, as those de- 
scribed in detail by Martyna et al.- A and Tuckermani 3 ^ 

Our calculations have been performed within the adia- 
batic (Born-Oppenheimer) approximation, which allows 
one to define a potential energy surface for the nuclear 
motion. An important issue in this kind of simulations is 
the proper description of interatomic interactions, which 
should be as realistic as possible. Since using true den- 
sity functional (DF) or Hartree-Fock-type calculations 
requires computer resources that would restrict enor- 
mously the size of our simulation cell, we obtain the Born- 
Oppenheimer surface from a tight-binding (TB) effective 
Hamiltonian, derived from DF calculations^ The TB 
energy consists of two parts, the first one is the sum of en- 
ergies of occupied one-electron states, and the second one 
is given by a pairwise repulsive interatomic potential. 36 
For the present study the H-H pair potential was tuned 
to reproduce the main features of known effective inter- 
atomic potentials, such as the Morse potential^ The ca- 
pability of TB methods to simulate different properties 
of solids and molecules has been reviewed by Goringe et 
aL— The convergence of the total energy with the sam- 
pling in reciprocal space was checked by using several 
sets of special fc-points^ We found that a set of 4 k- 
points provides already good convergence (relative error 
less than 0.001 % in the total energy). The use of only 
the r point introduces a small systematic error in the 
total energy that affects slightly the value of energy dif- 



ferences between different spacial configurations of H2 in 
silicon, with typical errors of about 0.01 eV. These results 
justify that the simulations presented in this work were 
performed by using only the T point for the reciprocal 
space sampling. 

Simulations were carried out on a 2 x 2 x 2 supercell of 
the silicon face-centered cubic cell with periodic bound- 
ary conditions, containing 64 Si atoms and a hydrogen 
(or deuterium) molecule. For comparison, we also carried 
out simulations of pure silicon, using the same supercell 
size. Sampling of the configuration space has been car- 
ried out at temperatures between 300 and 900 K. The 
electronic structure calculations were performed without 
considering a temperature-dependent Fermi filling of the 
electronic states, which is reasonable for this tempera- 
ture range. For a given temperature, a typical simulation 
run consisted of 10 4 PIMD steps for system equilibration, 
followed by 5 x 10 5 steps for the calculation of ensemble 
average properties. To keep a nearly constant precision 
in the path integral results at different temperatures, we 
have employed a Trotter number that scales as the in- 
verse temperature. In particular, we have taken PT = 
18000 K, which means P = 60 for T = 300 K. Quan- 
tum exchange effects between protons or deuterons were 
not considered, as they are negligible at the temperatures 
considered here, so that both atomic nuclei in a molecule 
were treated as if they were distinguishable particles. 

The simulations were carried out by employing a stag- 
ing transformation for the bead coordinates. The canon- 
ical ensemble was generated by coupling chains of four 
Nose-Hoover thermostats (with mass Q = (3h 2 /5P) to 
each degree of freedom^ To integrate the equations of 
motion, we used a reversible reference-system propagator 
algorithm (RESPA), which allows one to define different 
time steps for the integration of fast and slow degrees 
of freedom^ The time step At associated to the calcu- 
lation of DF-TB forces was taken in the range between 
0.1 and 0.4 fs, which was found to be appropriate for 
the interactions, atomic masses, and temperatures under 
consideration. For the evolution of the fast dynamical 
variables, including the thermostats and harmonic bead 
interactions, we used a smaller time step St = At/4. We 
note that for H2 in silicon at 300 K, a simulation run 
consisting of 5 x 10 5 PIMD steps needs the calculation of 
forces and energy with the TB code for 3 x 10 7 configura- 
tions, which has required the use of parallel computers. 



B. Calculation of anharmonic vibrational 
frequencies 

Vibrational frequencies of impurities in solids are im- 
portant characteristics, which depend on the site that 
they actually occupy and on its interactions with the 
nearby hosts atoms. In this context, the question arises 
whether the oscillator frequencies associated to an impu- 
rity can be extracted by assuming the host atoms fixed 
in the relaxed geometry corresponding to the minimum- 
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energy configuration. This is a method usually employed 
to calculate vibrational frequencies of impurities in crys- 
tals. On the other side, when the host atoms are allowed 
to relax by following the impurity motion, the potential 
energy surface is flatter than when the host atoms are 
fixed. To obtain an approach for the actual vibrational 
frequencies of the impurities, one can calculate the eigen- 
values of the dynamical matrix of the whole simulation 
cell, and obtain the frequencies in the harmonic approx- 
imation (HA). However, for light impurities the anhar- 
monicity can be appreciable, and the harmonic frequen- 
cies are only a first (maybe crude) approximation. 

To calculate anharmonic frequencies we will use here a 
method based on the linear response (LR) of the system 
to vanishingly small forces applied on the atomic nuclei. 
To this end, we consider a LR function, the static isother- 
mal susceptibility \ T > that is readily derived from PIMD 
simulations of the equilibrium solid, without dealing ex- 
plicitly with any external forces in the simulation. This 
approach represents a significant improvement as com- 
pared to a standard harmonic approximation. 3 — The ten- 
sor x T allows one to derive a LR approximation to the 
low-lying excitation energies of the vibrational system, 
that is applicable even to highly anharmonic situations. 
For a system with 3N vibrational degrees of freedom, the 
LR approximation for the frequencies reads 




where S n (n = l,...,3iV) are eigenvalues of x T , and 
the LR approximation to the low-lying excitation energy 
of vibrational mode n is given by Hui n . Details on the 
method and illustrations of its ability for predicting vi- 
brational frequencies of solids and molecules are given 
elsewhere^L^ 



III. RESULTS 
A. Minimum-energy configuration 

We first present results for classical calculations at zero 
temperature, where the atoms are treated as point-like 
particles without spatial derealization. The employed 
interatomic potential gives reliable results for molecular 
hydrogen in vacuo. The lowest-energy molecular config- 
uration corresponds to a distance Rq between hydrogen 
atoms of 0.74f A. At this distance we obtain for H2 in a 
harmonic approximation a stretching frequency of 4397 
cm -1 . 

For H2 as an impurity in silicon, we find a lowest- 
energy position for the center-of-gravity of the molecule 
located at an interstitial T site. The minimum energy is 
found for the H-H axis along a (100) crystal direction, 
with a distance between H atoms of 0.752 A. Moreover, 
changes in the energy for molecule rotation keeping its 
center-of-gravity at a T site are very small, in agreement 



with earlier calculations based on DF theor y 10 ! 11 ! 14 In 
silicon an increase in the H-H distance of about 0.01 A 
was found with respect to the molecule in vacuo, as ex- 
pected for an attractive interaction between each H and 
the nearby Si atoms. 

Assuming the H2 molecule at a T site, and oriented 
along the (100) direction, we find in the harmonic ap- 
proximation a stretching frequency of 4071 cm -1 , close 
to the harmonic value of 4015 cm -1 derived from the 
(anharmonic) vibrational frequency observed in Raman 
spectra 4 ^ This represents a decrease of more than 300 
cm -1 vs the harmonic frequency for the molecule in the 
gas phase, in line with a weakening of the H-H bond 
due to interaction with the silicon lattice, as discussed 
earlier. 10 ' 11 In the HA we also find frequencies uin — 954 
cm -1 and u>± = 1385 cm -1 (twofold degenerate) for mo- 
tion of the molecule along and perpendicular to the H-H 
axis in the silicon cage. These two vibrational frequencies 
are not expected to be observable because they will be 
mixed by the free rotation of the molecule (see below) . 



B. Kinetic energy 

We now turn to our PIMD simulations at finite temper- 
atures. To obtain insight into the motion of H2 around 
the tetrahcdral site, we will consider various models in 
which the number of degrees of freedom will be succes- 
sively increased. In particular, we will consider: (1) mo- 
tion of the H2 molecule in one dimension (along the H-H 
bond) in a fixed and unrelaxed silicon lattice; (2) free 
motion (in 3d) of the hydrogen molecule with fixed host 
atoms, and (3) free motion of H2 with mobile Si atoms. 
In the latter case, all 66 atomic nuclei in the simulation 
cell are treated as quantum particles. 

In our finite-temperature simulations for cases (2) and 
(3) , where the molecule can rotate around the T site, we 
observe a free molecular rotation, without any preferen- 
tial orientation. This is in agreement with earlier conclu- 
sions derived from theoretica l 10 ' 11 and experimenta l 12 ' 45 
works, and with the fact that the potential-energy surface 
for the rotation does not display deep minima. 

In Fig. 1 we show the kinetic energy of the hydrogen 
molecule in the three considered approaches. Symbols 
indicate results derived from our PIMD simulations using 
the so-called virial estimato r 40 ' 46 and solid lines represent 
the kinetic energy expected in a harmonic approximation. 
For Id motion of H2 (approach 1, squares) we find a slight 
increase in as temperature is raised. In this approach, 
results of the simulations are somewhat lower than those 
derived for the HA, as expected for a softening of the 
vibrations due to the anharmonicity of the interatomic 
potential. In fact, the linear-response method introduced 
in Sect. II. B gives in this case for the stretching frequency 
u) s = 3770 cm -1 at 300 K, which means a decrease of 
about 300 cm -1 with respect to the harmonic model for 
H2 in silicon (ui s = 4071 cm -1 ). 

Circles in Fig. 1 correspond to our approach 2 with H2 
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FIG. 1: Temperature dependence of the kinetic energy of 
molecular hydrogen in silicon for various approximations. 
Squares: motion of H 2 in one dimension with fixed host atoms; 
circles: free motion of H2 in a fixed silicon lattice; triangles: 
free motion of H2 with unrestricted motion of the Si atoms. 
Solid lines correspond to harmonic approximations for H2 mo- 
tion in one and three dimensions. Dashed lines are guides to 
the eye. 



FIG. 2: Temperature dependence of the vibrational part of 
the kinetic energy of H2 and D2, as derived from approach 3 
with free motion of all atoms in the simulation cell. Symbols 
indicate results derived from PIMD simulations: squares for 
H2 and circles for D2. Error bars are on the order of the 
symbol size. Dashed lines are guides to the eye. The dotted 
line corresponds to the classical limit with four degrees of 
freedom. 



moving in a fixed silicon lattice. Now we are dealing with 
six degrees of freedom, two of which correspond to molec- 
ular rotation. This approximation gives again values of 
the kinetic energy smaller than those predicted by the 
HA (solid line) . This harmonic approximation includes a 
classical description of the two rotational degrees of free- 
dom of the H2 molecule. To analyze the kinetic energy 
of the defect complex in model 3 (all atoms are free to 
move), we calculate for the simulation cell with and 
without the H 2 molecule: E k (defect) = E k (6A Si + H 2 ) 
- Ek (64 Si) . The results (triangles) lie appreciably above 
those obtained for model 2, indicating a nonnegligible 
coupling in the motion of interstitial molecule and host 
atoms. 

At the temperatures considered here, rotation of H2 
can be considered with a high precision to be classical. 
This means that its contribution to the kinetic energy 
of the molecule will be fc^T (two degrees of freedom). 
Then, we can subtract this classical energy from E k de- 
rived from the PIMD simulations to obtain a vibrational 
contribution to the kinetic energy E^ . This part of the 
kinetic energy is shown in Fig. 2 for H 2 (squares) and 
D 2 (circles). At low temperature it converges to val- 
ues close to 0.24 and 0.18 eV, respectively. This gives a 
ratio £7£(H 2 ) / £^(D 2 ) = 1.33, somewhat smaller than 
the limit 1.41 expected for harmonic vibrations at low 
temperatures. This deviation may be due to both anhar- 
monicity in the interatomic interaction and changes in 
the effective mass caused by coupling to the host atoms. 
This ratio decreases as T is raised, and amounts to 1.19 



at 900 K. For comparison we also present in Fig. 2 the 
kinetic energy corresponding to the classical limit with 
four vibrational degrees of freedom {2ksT, dotted line). 



C. Atomic derealization 

To study the spatial derealization of a quantum par- 
ticle from PIMD simulations, it is convenient to consider 
the center-of-gravity (centroid) of the quantum paths of 
the particle, defined as 

f =^X>> (2) 

i=i 

r, being the coordinates of the "beads" in the associated 
ring polymer. 

The mean-square displacement of a quantum particle 
along a PIMD simulation run is then given by: 

A 2 = ^X>-<r)) 2 ^, (3) 

where (...) indicates a thermal average at temperature T. 
After some straightforward manipulations, one can write 
A 2 as 

A^A^ + A 2 ,, (4) 
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FIG. 3: Spatial derealization of atomic nuclei (protons) in 
H2. Diamonds indicate the mean-square displacement of the 
centroid of the quantum paths, Aq, and squares correspond 
to the mean-square radius-of-gyration of the paths, Aq. 

with 
and 

A 2 c = ((r-(r» 2 ) . (6) 

The first term, Aq, is the mean-square "radius-of- 
gyration" of the ring polymers associated to the quantum 
particle (atomic nucleus) under consideration, 1 - This is 
a measure of the average extension of the paths and, 
therefore, of the importance of quantum effects in a given 
problem. The second term in Eq. (j4]) is the mean-square 
displacement of the center of gravity of the paths. This 
term is the only one surviving at high temperatures, since 
in the classical limit each path collapses onto a single 
point (hence with a vanishing radius-of-gyration). For 
situations in which the anharmonicity is not extremely 
large, the distribution of r is similar to that of a classi- 
cal particle in the same potential, and thus Aq can be 
considered as a kind of semiclassical derealization. 

Going back to our problem of H2 in silicon, for each 
hydrogen atom in the molecule we have calculated sep- 
arately both terms giving the atomic derealization in 
Eq. (j4|). Shown in Fig. 3 are the values of Aq (spreading 
of the quantum paths, squares) and A c (centroid der- 
ealization, diamonds), as derived from our PIMD simula- 
tions at several temperatures. In this plot, one observes 
that Aq is much larger than Aq in the whole tempera- 
ture range under consideration. This is not strange if one 
takes into account that the molecular rotation around the 
interstitial T site can be considered as a classical motion 



at these temperatures. In fact, the order of magnitude of 
this spatial derealization can be obtained from the free 
motion of a particle on a spherical surface with radius 
equal to half distance between atoms in an H2 molecule. 
For a distance R = 0.77 A, we obtain a mean-square clas- 
sical displacement of 0.15 A 2 , close to the value of A^ 
at 300 K. This magnitude increases for rising tempera- 
ture, as expected for an increase in the fluctuations of 
the distance from each H atom to the average position 
(T site). 

For the spreading of the quantum paths of each H atom 
we obtain at room temperature Aq = 0.03 A 2 , and it de- 
creases as temperature is raised. This gives for the paths 
an average extension of ~ 0.1 A at 300 K, much smaller 
than the H-H distance, and thus justifying the neglect 
of quantum exchange between protons. Moreover, the 
fact that Aq is much smaller than A c in the tempera- 
ture range considered here does not mean that quantum 
effects arc irrelevant, but is a consequence of the enhance- 
ment of A% due to molecular rotation. 



D. Interatomic distance 

As mentioned above, the interatomic distance between 
hydrogen atoms increases when the molecule is intro- 
duced from the gas phase into a silicon crystal, due 
to the interaction between H and host atoms. For the 
minimum-energy distance we found Rq — 0.752 A, which 
is smaller that the values obtained in earlier calculations 
(0.788 A in Ref. and 0.817 A in Refs. [ll]|47|). 

We now present the temperature dependence of the 
mean distance H-H for the three approaches considered 
above to study molecular hydrogen in silicon. Our results 
are displayed in Fig. 4, where symbols represent data 
points derived from PIMD simulations. For approach 1 
(Id motion in a fixed lattice), we find at 300 K a mean 
distance R — 0.780 A, which represents an appreciable 
increase vs the distance obtained for the minimum-energy 
configuration. In this model, R increases very slowly as a 
function of T (squares in Fig. 4), since molecular rotation 
is not allowed and the molecule expansion is only due to 
the increasing population of excited vibrational levels. In 
fact, we find dR/dT = 1.0 x 10~ 6 A/K. 

When molecular rotation is allowed in a fixed lattice 
(approach 2), we observe an increase in R (see circles in 
Fig. 4). At 300 K we found R = 0.785 A, about 5x 10~ 3 A 
larger than for Id motion. Now R rises with temperature 
much faster than in approach 1, with dR/dT = 1.1 x 10 -5 
A/K. Next, we allow the Si atoms to move, introducing 
the full quantization of all degrees of freedom in the sim- 
ulation cell, and we obtain a reduction of the distance H- 
H with respect to approach 2. This can be understood 
as due to a softening of the effective Si-H interaction, 
which decreases as a consequence of Si motion. In this 
case with full motion of the 66 atoms in the cell, we find 
dR/dT = 1.6 x 10" 5 A/K, which means a larger slope 
than in approach 2. 
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FIG. 4: Mean distance between H atoms in an H2 molecule 
in silicon. Symbols correspond to different approximations for 
the molecular motion. Squares: motion of H2 in one dimen- 
sion with fixed host atoms (approach 1); circles: free motion 
of H2 in a fixed silicon lattice (approach 2); triangles: free 
motion of H2 and host atoms (approach 3). Error bars are in 
the order of the symbol size. Dashed lines are linear fits to 
the data points. 
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FIG. 5: Mean interatomic distance for H2 and D2 molecules 
in silicon, as a function of temperature. Symbols indicate re- 
sults derived from PIMD simulations for approach 3, in which 
all atoms are mobile: squares for H2 and circles for D2. Error 
bars are on the order of the symbol size. Dashed lines are 
linear fits to the data points. 



E. Stretching frequency 



It is interesting to compare these changes in the mean 
distance R with those corresponding to molecular hy- 
drogen in the gas phase. To this end we have car- 
ried out some PIMD simulations of an isolated hydrogen 
molecule with the same interatomic potential. In this 
case we obtain an increase in R with temperature given 
by dR/dT = 7. 5xlQ~ 6 k/K, a value clearly smaller than 
those obtained for H2 in silicon in our approaches 2 and 
3. This means that, for H2 in silicon, the change of inter- 
atomic distance with temperature is controlled by both 
the centrifugal expansion due to rotation, and interaction 
with the nearby host atoms. 

PIMD simulations can be also employed to study the 
isotopic dependence of the mean interatomic distance 
R. The molecular expansion with respect to the lowest- 
energy classical geometry is due to a combination of 
anharmonicity with quantum delocalization. One ex- 
pects smaller distances for molecular deuterium due to 
its smaller vibrational amplitudes. In fact, at 300 K we 
found for D2 in silicon, R = 0.767 A, to be compared 
with R = 0.776 A for H2 at the same temperature, and a 
distance Rq = 0.752 A for the lowest-energy position in 
the classical limit. In Fig. 5 we present the temperature 
dependence of the mean distance for both H2 and D2 , as 
derived from our PIMD simulations for approach 3 (full 
motion of molecular hydrogen and host atoms). For D2 
we find dR/dT = 1.5(1) x 10~ 5 A/if, which coincides 
within error bar with the slope obtained for H2 in silicon 
in the temperature region from 300 to 900 K. 



The stretching frequency of H2 is an important finger- 
print of the molecule, that in fact has been used to de- 
tect and characterize this impurity in the silicon bulk^^ 
This stretching vibration has been found at 3618 cm -1 
(at 4 K) independently by Raman 8 and infrared absorp- 
tion spectroscopies^ 

In Fig. 6 we show the temperature dependence of u) s 
for H2 in the three approaches considered here, as derived 
from the LR method presented above. In approach 1 (Id 
motion) the frequency decreases slightly in the analyzed 
temperature range. In approaches 2 and 3, the coupling 
between molecular rotation and vibration causes an ap- 
preciable change of uj s with the temperature. For model 
2 (fixed Si lattice) we find dw s /dT = -0.13 cm -1 / K > to 
be compared with a slope of dcu s /dT — —0.24 cm. _1 /K 
for model 3, which includes motion of the host atoms. 
Thus, motion of the Si atoms causes a significant change 
in duj s /dT, which becomes almost twice larger than in the 
case of a static Si lattice. It is interesting that at room 
temperature u> s is smaller for model 2 than for approach 
3, but due to its faster decrease in the latter approach, 
uj s becomes smaller for model 3 at high T. 

Something similar has been obtained for the stretch- 
ing vibration of D2. In particular, for approach 3 we find 
a rather constant ratio between the stretching frequen- 
cies of H2 and D2, that amounts to 1.37(1), somewhat 
smaller than the ratio expected in a harmonic approxi- 
mation (1.41). Experimentally, a ratio of 1.37 is found 
from infrared^ and Raman 8 -^ spectra of H2 and D2 in 
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FIG. 6: Frequency of the stretching vibration of the H2 
molecule in silicon as a function of temperature. Symbols rep- 
resent results derived from PIMD simulations in three differ- 
ent approximations: squares, motion of H2 in one dimension 
in a fixed silicon lattice; circles, free motion of H2 with fixed Si 
atoms; triangles, free motion of H2 and Si atoms. Error bars 
correspond to the statistical uncertainty in the PIMD simu- 
lations. A black diamond indicates the stretching frequency 
obtained by Raman spectroscopy at room temperature.— 



silicon, a little smaller than the ratio 1.39 observed for 
these molecules in the gas phased 

For the HD molecule in silicon, an infrared study al- 
lowed to determine the energy of the first excited rota- 
tional level 4^ In fact, a value of 73.9 cm -1 was found 
for the wave-number difference between the levels J = 
and J = 1, somewhat lower than that corresponding to 
the gas phase (89.3 cm -1 ). By scaling that wave-number 
difference with the reduced mass, we expect for H2 an 
energy difference of about 99 cm -1 . Since our PIMD 
simulations yield results for the average frequency u> s , 
one can estimate a frequency shift from the rotational 
energy, taking into account the population and degen- 
eracy of the different levels^ By considering only the 
levels J = and J = 1, one would expect at room tem- 
perature a frequency shift dw s /dT on the order of —0.05 
cm /K, clearly lower than the value found from our 
simulations for approach 3 (—0.24 cm /K). This is not 
strange, taking into account that at these temperatures 
higher rotational levels will be excited, contributing to a 
larger decrease in the average frequency. However, the 
actual position of these levels further than J = 1 is not 
known at present, and a more detailed comparison with 
our results is not possible. 

We note that the quantum treatment of atomic nuclei 
in molecular dynamics simulations is crucial to give a re- 
liable description of the vibrational frequencies of light 
atoms like hydrogen. In fact, we have applied the LR 
method to calculate the stretching frequency lu s from 



classical simulations. At 300 K we found for H2 in silicon 
a frequency uj s = 4039 cm -1 (for full motion of inter- 
stitial hydrogen and host atoms), to be compared with 
u) s = 3728 cm -1 derived from PIMD simulations. As ex- 
pected, the classical value is much closer to the frequency 
u> s = 4071 cm -1 obtained in a HA for H2 in silicon. 



IV. DISCUSSION 

In Sec. Ill we have presented results of our PIMD sim- 
ulations for H2 and D2 in silicon. The main advantage 
of this kind of simulations is the possibility of calculat- 
ing energies at finite temperatures, with the inclusion of 
quantization of host-atom motions, which are not easy 
to be accounted for in fixed-lattice calculations. Isotope 
effects can be readily explored, since the impurity mass 
appears as a parameter in the calculations. This includes 
the consideration of zero-point motion, which together 
with anharmonicity gives rise to non-trivial effects. In 
addition, the vibrational motion of H2 is coupled with 
molecular rotation, leading to a change in the stretching 
frequency with temperature. 

As mentioned above, an important feature of isolated 
H2 molecules in semiconductors is their stretching vibra- 
tion u> s . In a harmonic approximation, the tight-binding 
potential employed here yields for H2 in silicon a fre- 
quency uj s = 4071 cm" 1 vs 4397 cm -1 for H2 in the gas 
phase, which means a reduction of about 330 cm -1 due 
to interaction with the host atoms. This reduction is ac- 
companied by an increase in the H-H distance, as shown 
in Sect. III. C. An additional decrease in uj s is obtained 
when anharmonic effects are taken into account in a one- 
dimensional motion of the molecule in a fixed lattice. In 
fact, at 300 K the LR calculations give in this case co s 
— 3770 cm -1 , which means a decrease in frequency of 
about 300 cm" 1 with respect to the HA for H2 in silicon. 
This frequency change due to anharmonicity is in the 
order of that derived in Ref. |47J from DF calculations, 
namely Aw s = —408 cm . This frequency is further 
lowered when full (quantum) motion of the molecule and 
host atoms are allowed, giving uj s = 3728 cm . In this 
latter reduction there is a contribution of two competing 
effects: coupling between molecular rotation and vibra- 
tion, and interaction with Si atoms, whose motion allows 
for a larger delocalization of the H2 molecule in the in- 
terstitial space. 

In general we observe a correlation between oj s and 
mean interatomic distance R in the H2 molecule, in the 
sense that a rise in oj s is accompanied by a decrease in 
R. This is in line with the general trend found by Van 
de Walleii for molecular hydrogen in crystalline semi- 
conductors, as derived from DF calculations at T = 0. 
However, this trend is not so strict when atomic motion 
is included at finite temperatures, as derived from Figs. 4 
and 6. In this case, the mobile Si atoms may contribute 
to an additional decrease in the stretching frequency of 
H2 by a rise in the effective mass associated to this vi- 
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brational mode. 

In connection with this, it is clear that theoretical tech- 
niques to deal with the electronic structure of solids have 
been improving their precision over the years. For var- 
ious purposes, the accuracy currently achieved by these 
methods is excellent, when comparing their predictions 
with experimental data. However, quantum nuclear ef- 
fects limit the accuracy of state-of-the-art techniques to 
predict actual properties of light impurities in solids. The 
answer to this question has to be found in ab-initio path- 
integral simulations, where both electrons and nuclei are 
treated directly from first principles. But even in this 
case the question is not simple when one has to deal with 
phenomena such as molecular rotation at low tempera- 
tures, where a proper description of quantum rotation 
has to be included in the formalism. 

There is an important challenging point that should be 
considered in future work. It refers to considering cou- 
pling between nuclear spins in the hydrogen molecule, i.e. 
dealing separately with ortho and para-H2 (both have 
been observed in silico n 12 i 44 i 51 ). This becomes specially 
relevant at low temperatures, where the quantum charac- 
ter of molecular rotation has to be explicitly considered 
in the simulations. Usually this kind of calculations have 
been carried out by assuming the molecule to be a rigid 
rotor, without taking into account vibrations and defor- 
mations, and thus neglecting the ro- vibrational coupling. 



An analysis of hydrogen diffusion in silicon is out of 
the scope of this paper. Actual diffusion coefficients are 
not directly accessible with the kind of simulations em- 
ployed here, since the time scale employed in the calcula- 
tions is not readily connected to the real one. In this 
respect, PIMD simulations could be applied to study 
quantum diffusion of H2 in silicon, by calculating free- 
energy barriers in a way similar to that employed earlier 
to study the diffusion of atomic hydrogen in metals 2 ^ and 
semiconductors j 24 ' 52 

In summary, the PIMD method has turned out to be 
well-suited to study finite-temperature equilibrium prop- 
erties of hydrogen molecules in silicon. This has allowed 
us to notice the importance of anharmonicity and ro- 
vibrational coupling in order to give a realistic descrip- 
tion of the properties of these interstitial impurities. An- 
harmonicity shows up in the stretching motion of the 
molecules, causing important shifts with respect to the 
harmonic expectancy. 
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